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ABSTRACT 



To help understand the stability of cold, viscous boundary layers in 
geophysical contexts such as lava lakes and mantle convection, the following 
model problem is analyzed: Beneath a shear-free horizontal boundary, a thin 

layer of very viscous fluid overlies a deep layer of less viscous, less dense 
fluid. The initial unstable equilibrium is perturbed, and the growth of the 
disturbance is followed, including the nonlinear effects of large amplitude, 
by a long^wave analysis. The result shows that in the final catastrophic 
growth the peak thickness of the upper layer approaches infinity inversely 
proportional to the remaining time. (This conclusion also applies to fluids 
with power-law rheology.) Thus nonlinear effects greatly enhance growth. 
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Summary of notation 



Pi, Ui, P2, P2, A P> g, x > z > t, u, w : see Figure 1 

P : local density ( Pi or P 2 ) 

v : either kinematic viscosity 

k : characteristic wavenumber 

di, 62 : equilibrium depths of each layer 

6(x,y,t) : interface location 

P : local pressure (Pi or P 2 ) 

p : local reduced pressure (reduced by pg(z-di)) 
c* : stress tensor 

T ij : deviatoric stress tensor (without pressure contribution) 

a - 2 / U 1 : viscosity ratio 

3 = d 2 / d 1 : depth ratio 

e = kdi : dimensionless wavenumber 

Fj_j : reduced force tensor, 2D 

F(t) : reduced force for ID disturbance 

d e (t) : current "equilibrium” thickness 

xq : initial position of fluid cross section 

<5o(xo) : initial thickness profile 

a(x,t) : disturbance amplitude (dimensionless, = 6-1) 
a o(xo) : initial amplitude 

t* : singular time for any cross section to reach infinite thickness 
t = t-t* : time relative to singular time 

L : length of layer, or wavelength of periodic disturbance 
Lq : initial length or wavelength 

n(x 0 ,t), f (n) , 5(n) : f or similarity solutions describing plume 

x*(t) : "initial” position of fluid section currently at plume position 



a : arbitrary power for peak profile 
e : small remainder 
0 : linearized growth rate 

t(6) = t*-t : time until singularity for power-law fluid layer 

a(t), A(t), 5(x,t), f(5), g(0 : for large-amplitude similarity solution 
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1. INTRODUCTION 

This work describes the circumstances under which the motion of a 
buoyantly unstable horizontal film of viscous fluid is limited by normal 
viscous stresses, and gives an analysis showing how the growth of disturbances 
to the interface becomes greatly enhanced when the disturbance amplitude 
becomes large, leading to the formation of plumes in a finite time. A 
companion work will show that in certain cases this same dynamic balance 
controls the structure of the cold boundary layer in Benard convection with 
strongly temperature-dependent viscosity. The behavior described here thus is 
of geophysical relevance to such flows as those in lava lakes, the thermal 
convection of planetary mantles, and possibly also in convection in the 
Earth’s solid core. The Rayleigh-Taylor problem considered here is the 
simplest example of this dynamic balance. 

The problem is illustrated in figure 1. Two layers of immiscible fluids 
are confined between shear-free boundaries, with the upper fluid layer being 
denser and much more viscous than the lower, and both fluids are so viscous 
that neither inertia nor surface tension is significant. A linearized 
analysis (appendix A) shows that the most unstable wavelength is long compared 
to the upper layer thickness. And scaling the problem shows that this fastest 
growth occurs in a broad range of wavelengths for which the lower fluid is 
effectively passive and hydrostatic. Thus we give a long-wave analysis, where 
the lower fluid is treated as inviscid, that exploits the fact that, for long 
waves, the disturbance amplitude can get very large (compared to the average 
layer thickness) while the slope of the interface remains small; this follows 



from mass conservation. 
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The results of this analysis show the nonlinear effects of large 
amplitude. In particular, it predicts that thickness maxima in the layer will 
reach infinite thickness in finite time, with the final catastrophic growth 
inversely proportional to the time remaining before the singular time. This 
conclusion is not unique to a Newtonian rheology; the corresponding analysis 
for a power-law fluid (appendix B) shows the same catastrophic inverse-time 
growth of large peaks. 

The small-slope analysis predicts that at the singular time, a peak will 
locally have the shape ® 01 |x|“ 2/3 , but that afterward, as the plume acts as a 
sink of fluid, the shape will change to ^ “ |x|~ 1/2 . 

Of course the small-slope approximation must fail before this, but we 
argue that the failure is only local; the slope remains small everywhere 
except in an asymptotically small neighborhood around each developing 
singularity. Furthermore, these small regions where the small-slope 
approximation breaks down have only an asymptotically small effect on the 
dynamics of the rest of the layer. Therefore the small-slope equations will 
continue to apply almost everywhere, even after the singular time, when peaks 
form downwelling plumes. (There is a family of similarity solutions, given in 
appendix C, that describe plume behavior.) Thus the same equations describe 
the disturbance from the initial linear growth through ’the rapid large- 
amplitude growth to the final draining of the layer by the plumes. 

In the companion paper on the analogous thermally driven buoyant 
instability, there is a family of steady solutions (with plumes), and we use 
the small-slope equations to follow the development from initial conditions 
through to the final steady state. 
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2. PROBLEM STATEMENT 

Two horizontal layers of distinct Newtonian incompressible fluids are 
bounded above and below by horizontal shear-free boundaries (see figure 1). 

The upper fluid (of density Pi and viscosity V\) is denser and much more 
viscous than the lower fluid (of density P2 < Pi and viscosity ^2 << Ui), and 
the upper layer is very thin (di << d2). Both fluids are assumed to be so 
viscous that we can neglect both surface tension and inertia. The latter 
assumption requires 

Ap gL 3 /PV 2 « i 

where L is the largest length scale, ^p = P 1 -P 2 , g is gravity, P is either 
density, and v is either kinematic viscosity; this is easily satisfied in 
mantle flow. 

Initially both fluids are at rest in unstable equilibrium. Then at t = 0 
the interface is slightly disturbed; thereafter the interface position is 
given by ^(x,y,t). (Later we will assume that the disturbance varies only in 
the x direction; a one-dimensional disturbance allows some remarkable 
simplifications, and illustrates the physical behavior more clearly.) 

A reduced pressure p is defined by: 

P( x , y , z , t ) = p (x , y , z , t ) + Pg(z-di) (2.1) 

where p is the local density ( p 2 or Pi). Then the governing equations are: 

v p = U v2 u (2.2a) 

v 'u = 0 (2.2b) 



and the boundary conditions are: 
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= 0 and at z = (dl+d2): w = u z = 0 


( 2 . 3a , b) 


= 5 (x,y, t) : U = 0_, 


(2.3c) 


[n°] = A Pg(6-dl)n, 


(2.3d) 


\ + u + v ^y = w 


(2.3e) 



where LJ indicates the change in the enclosed quantity across the interface 
(downward), is the reduced stress tensor, and the unit normal to the 
interface n has components: 

n = ( -fi x , -« y> 1 )/* / (l+« x 2 +«y 2 ) (2.4a) 

Physically, the conditions (2.3a-d) are that the boundaries are impermeable 
and exert no shear, and across the interface both velocity and stress are 
continuous, while (2.3e) is the kinematic condition that the interface moves 
as a material surface. 

These equations and boundary conditions, along with specification of the 
initial interface position ^(x,y,0), define the problem without approximation. 
To facilitate analysis, however, we limit attention to the development while 
the slope of the interface remains small, i.e., and are both small. 

Then in the interface stress conditions we can approximate the unit normal, to 
lowest order, by: 

n = ( -fi x , -«y, 1 ) (2.4b) 

and neglect any resulting terms quadratic in and ^y.' It should be noted 
that the resulting small-slope equations still retain the leading nonlinear 
terms due to slope and are still applied at z = <$(x,y,t), in contrast to the 
linearized conditions used for the small-amplitude analysis (appendix A). 

The small-slope approximation applies trivially to small-amplitude 



disturbances. But more importantly, the small-slope approximation applies to 
long-wavelength disturbances even when the amplitude becomes large (until the 
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amplitude is comparable to the wavelength). Thus we can use the simplified 
equations to analyze the large-amplitude, nonlinear effects in long-wave 
instabilities, which we do in the next section. 

3. LONG-WAVE ANALYSIS 

The disturbance can grow in a variety of different ways, depending on the 
viscosity ratio, the depth ratio, and the dimensionless wavenumber: 

a = U 2 /Ui, 6 = d 2 /d i , e = kdi (3.1) 

(where the disturbance has a characteristic wavenumber k). The linearized 
analysis in Appendix A gives the small-amplitude growth rates for the full 
range of all parameters. A scaling analysis (see Canright, 1987, App. B) 
shows that the same growth regimes apply even for large amplitudes, as long as 
the slope of the interface remains small. 

Here we focus on the case where the lower fluid layer is much less 
viscous and deeper than the upper, so: 

° « 1, 8 » 1 (3.2) 

Then there is a range of wavelengths, which includes the fastest growing 
wavelength, over which the growth rate is nearly constant (see fig 6): 

max( a , *'( a /8) ) « E « 1 (3.3) 

In this range, the growth is limited by normal stresses in the upper fluid, 
which moves nearly horizontally, while the lower fluid 'is passively moved by 
the interface. Outside this range, for waves short compared to dl , the growth 
is reduced because only a fraction of layer 1 is mobilized, and for long 
enough waves the viscous resistance of fluid 2 slows the growth. 

We examine the finite-amplitude growth of long-wave disturbances in this 
fas test-growth regime, exploiting the fact the slope of the interface remains 
small even for large amplitudes (until the disturbance grows to the order of 
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the initial wavelength). What follows is the leading-order asymptotic 

4 

analysis; since the growth rate is constant to within 0( e ), we cannot expect 
this analysis to predict the single wavelength giving maximum growth. 

For wavelengths in this range, the interface motion is controlled by the 
dynamics of the upper fluid, while the lower fluid is effectively inviscid and 
hydrostatic; the error we make in neglecting the viscous resistance of the 
lower fluid is much smaller than 0( e ). Hereafter we refer only to fluid layer 
1 and drop the subscript 1 where clear. The motion is quasi-horizontal as 
both surfaces of the layer see no shear. It is driven by the horizontal 
variations of the buoyant pressure and resisted by the normal viscous 
stresses . 

We now estimate the scales of the different terms, assuming for this 

purpose that the disturbance varies only in the x direction. Then there are 

two length scales, x ~ k" 1 , z ~ dl, and ~ e . From continuity (and the 

impermeable boundary), w/u ~ e . Consequently, comparing vertical and 

horizontal momentum equations shows that p 2 /p x ~ and the interface shear 

stress condition (2.3c) implies u z ~ w x , so u 2 /u x ~ e . The last scaling in 

2 

turn implies that the variation of u across the layer is 0( e ) smaller than u 
(similarly for p), though u 2Z ~ u xx . 

When the initial disturbance to the interface varies in two dimensions, 
over a length scale long compared to the thickness of the upper layer but not 
so long that the resistance of the lower fluid becomes important, then the 
above scaling still applies. That is, u, v and p are independent of z, and: 

w = -z(u x +Vy) (3.4) 

all within an error of 0(e 2 ). From the normal stress condition, to 0(e 2 ): 



"P + T zz = Ap>g( 6-d 1 ) 



(3.5) 
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where is the deviatoric stress tensor. 

Consider a force balance in the x direction on a small column of layer 1 

(see figure 10). We reduce the force on all surfaces by the hydrostatic 

pressure gradient due to fluid 2; this does not affect the balance. Recalling 

(2.1) and from (3.5), the difference in total pressures is: 

P 1 -P 2 = p + Apg(z-di) = T zz - Apg(6-z) (3.6) 

Then the reduced force balance involves only the sides of the column, since 

the boundary is shear-free and the reduced force on the interface is zero, so: 
Ay A x 

— F xx ( x ) ] dy + J fF X y(y+Ay) _ F X y(y)] dx = 0 (3.7) 

y 6 x 

F xx ~ ^ [-(P1-P2) + T xx^ [ A pg 6 / 2 -h ^( T xx - T zz )l 

°6 

F X y “ / T X y dz [^ T xyl 

0 

where F^ is a 2-D tensor representing the reduced force in the layer: any 
vertical plane through the layer defines a horizontal normal direction, and 
the product of that normal and the tensor F^j gives the reduced force vector 
acting on that plane. Dividing by A x Ay an d taking the limit as A x ,Ay + 0: 

3/3x(F xx ) + 9 / a y(F xy ) = 0 (3.8) 

Indeed, a force balance in the y direction shows: 

=0 (3.9) 

F U - *<-P + T ij> 

P E - A Pg 6 /2 + T zz 

where P is the average (reduced) pressure and is the Kronecker delta. 

The kinematic interface condition becomes: 

D^/Dt = + u^ x + v^y = - ^ ( u x + Vy) (3.10) 

These two equations (3.9-10) govern the growth of long-wave disturbances. 
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Note that the above derivation is independent of the rheology of the fluid 
layer; any constitutive relation could be specified. 

The analysis simplifies considerably for a one-dimensional disturbance; 
we develop only this case. In general, the reduced force F^j is a two- 
dimensional tensor with zero divergence, reflecting the balance of buoyant and 
viscous forces in the layer. However, when the disturbance is one- 
dimensional, then only one component of F^j is nonzero; we call this scalar F. 
Then (3.9-10) become: 

9 / 9 x F = 0 (3.11) 

D 6 /Dt = 6 t + u<$ x = -<5u x (3.12) 

the former of which integrates directly to give: 

4<$u x + (Apg/y)6 2 /2 = F(t) (3.13) 

This says that the force on any cross section of the layer, reduced by 
the hydrostatic pressure force due to the lower fluid, is the same everywhere 
in the layer, i.e., F(t). (This result is not surprising, in that the only 
outside force on the layer is just the external hydrostatic head.) F(t) 
depends on the conditions specified at the ends of the layer, e.g., if the 
layer had abrupt ends with no applied force, surrounded by fluid 2, F would be 
zero. Arbitrary end conditions can thus be incorporated via F(t). 

We will treat three cases where F is simple to evailuate. In the first 
case, the layer is of infinite extent, but the disturbance is localized, so 
that far away the layer remains in (unstable) equilibrium, and F is constant. 
This case has a simple analytic solution that illustrates clearly the 
nonlinear effects of finite amplitude. The second case is a periodic 
disturbance (or equivalently, a layer with shear-free end walls). Here, the 
requirement that u = 0 at both ends of a period determines F(t) as an integral 
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property of the shape of the layer. In the last (and simplest) case the layer 
is of finite extent so F = 0; this solution also applies approximately when ^ 
is large and F can be neglected. 

In dimensionless form, the governing equations (3.12-13) become: 

D^/Dt = + u^/^x = -$9u/9x (3.14) 

S9u/3i + l 2 = F(t) = d e 2 (t) (3.15) 

which can be combined to eliminate u by adopting a Lagrangian formulation, 
where the fluid ’’particle” in this case is a material cross section of the 
layer, to give: 

D*/Dt = 3 2 - d e 2 (t) (3.16) 

Here D/Dt - <*/3t + u<V^x is the time derivative following a fluid cross 
section, ^ - <Vdl, t = °t , 0 = ^Pgdl/(81 J ) is half the small-amplitude growth 
rate, x = kx, u - ku/°, k is the wavenumber, and d e (t) = [ F( t ) / ( ^Pgd 1 2 /2 ) ] 1 7 2 
is the current (nondimens ional ) equilibrium thickness, i.e., that thickness 
with no tendency to grow or shrink, so d e (0) = 1. Equation (3.16) is a 
Ricatti equation, and can be integrated analytically for several different 
forms of d e (t), or numerically for arbitrary d e (t). 

This result shows that the vertical motion of the interface depends onlv 
on the local layer thickness and the current equilibrium thickness (which may 
depend on the overall shape of the whole layer). Thus any fluid cross section 
does not care what its immediate neighbors are doing, and the growth is 
insensitive to wavelength, as expected. Of course, the horizontal motion of 
any cross section depends on the shape of the whole layer. 
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4. SOLUTIONS 

4.1. LOCALIZED DISTURBANCE 

The most illustrative case is an infinite layer with a disturbance of 
finite extent. Then the force on the ends of the disturbed region remains 
constant: d e = 1. This has a simple solution (drop tildes): 

6<x 0>t ) = <1 + a(xq ) = (4.1) 

(1 - A(x 0 )e2t) (6o( x o^ + U 

where 6o(xg) = 6(xg,t=0) and x 0 identifies a particular fluid cross section by 
its initial position. This becomes singular in finite time: when 
t = t* = In | 1 / A | the thickness 6 goes to °° or 0, depending on whether 6 q was 
greater or less than 1. (When the minimum thickness in the layer first 
reaches zero, then d e becomes zero, and the above solution becomes invalid for 
the entire layer.) 

When the initial perturbation amplitude a 0 (x Q ) = <S Q (x 0 )-l is small, the 
solution reduces to: 

6 53 — — + 0 (ag 2 e 2 t) (4.2a) 

1 - a 0 e2t/ 2 

and while the amplitude a(x Q ,t) = 6(x 0 ,t)-l remains small: 

a - a 0 e2t Q + a 0 (e2t-i)/ 2 ) (4.2b) 

Initially this gives the exponential growth of the linearized solution, and as 
nonlinearities become important the perturbation growth accelerates for peaks 
(a>0) and retards for troughs (a<0). As a result of this and volume 
conservation, the peaks get narrower and sharper while the troughs broaden and 
flatten. We can roughly estimate the duration of the linear behavior by when 
this small-amplitude solution gives a ~ 1: At ~ In | 1 / a q | . 
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As the amplitude gets large, the growth becomes algebraic in the time 
remaining before the thickness becomes singular («> or 0): 

6 0 > 1: 6 » l/(t*-t) (4.3a) 

6 o < 1 : 6 « t*-t (4.3b) 

which shows that the rapid nonlinear growth occurs over a time At ~ 1 , or, 
dimensionally, 8 y/ Apgd 1 . 

The catastrophic growth shown by the inverse relation (4.3a) between peak 
thickness and remaining time is strikingly different from the exponential 
growth of small amplitudes. For large peaks the relevant time scale is the 
time remaining before the singularity, which is inversely proportional to the 



current dimensional peak thickness 6, 



max • 



p/Apg 6 



max • 



This shows that 



large-amplitude effects drastically enhance the growth of peaks. 

To see the shape of the layer requires Eulerian coordinates. As x is the 
integral of the strain dx/dx 0 , and because volume in the layer (and in each 



fluid cross section) is conserved, dx/dx 0 = 6q/6: 

*0 

x(x 0 ,t) = / 6 0 (x 0 )/<5(x 0 , t) dx 0 

0 

where the x origin is chosen at some stationary fluid cross section. 



(4.4) 



We find that, at least initially, the overall strain increases, i.e., the 
perturbed section stretches out in the x direction, for any initial 
perturbation with a zero mean. To see this, we take the time derivative of 



the position L of the end of the disturbed region: 

L 0 

DL/Dt — J 6 q( 1- 6 2 ) / 6 2 dxQ 
0 

Then at t = 0, in terms of the amplitude a(x,t): 

L L 

a 2 / 6 dx > 0 



(4.3a) 



DL/Dt = / (-2a + a 2 / 6) dx = ^ a 2 /6 dx > 0 (4.5b) 

since fa dx = 0. That the disturbed region stretches out in the x direction 
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is surprising, but is entirely consistent with the assumption that the force 
on the ends remains constant. 

For example, if the initial perturbation is sinusoidal we get: 

6 0 = 1 - b cos Xq (A. 6a) 



6(x 0 , t) 
x(x 0 ,t) 
+ 



2 - b (l+e 2t ) cos x n 
2 - b (l-e 2t ) cos x q 

= bC sin x 0 + (D-C)xq 

2CD _____ /(E 2 -b 2 ) tan(x n /2) 
arctan “ 



(4.6b) 



(4.6c) 



/(E 2 -b 2 ) E - b 

C = (e 2t -l )/ (e 2t +l ) , D = 4e 2t /(e 2t +l) 2 , E = 2/(e 2t +l) (4.6def) 

and the overall strain is D-C + CD//( E 2 -b 2 ) , which increases monotonically 
with time. This solution is graphed as 6(x) for various times in figure 2. 
When the initial amplitude b is small, the solution simplifies: 



6 « 1 - t A - c .£_ s -_ x -a , A ( t ) = b e2t/2 (4.6g) 

1 - A COS Xq 

x «= - tan“l[ — * — tan(x Q /2) ] - x Q (4.6h) 

/ ( 1-A 2 ) J ( 1+A) 

4.2. CONSTANT WAVELENGTH DISTURBANCE 

For the second case, we consider the more physically reasonable situation 

of a layer of finite extent bounded by fixed (shear-free) end walls (or the 

equivalent situation of a periodic disturbance of infinite extent). Then the 

end forces must vary with time to give u = 0 at both ends. Integrating (3.9) 

in x shows that fixed ends require: 

L L 

d e 2 (t) * (/ <5 dx) / (/ 1/6 dx) 

0 0 (4.7; 

where L is the entire length of the layer (or one wavelength), so the 



numerator is just the total volume of fluid in the layer, a constant. Thus 



13 



the thinnest parts of the layer have the greatest effect on d e (t). We find 

that, as a result of (A. 7), d e (t) decreases in general, although the effect is 

small for small amplitude. As the disturbance grows and the troughs deepen 

and broaden while the peaks become narrower, d e can only decrease with time. 

Combining (4.7) and (3.10) (using dx = (6 0 /6)dx Q ) gives a single partial 

integro-dif ferential equation that can be solved numerically for 6(x Q ,t) given 

an arbitrary initial profile 6q(x q ): 

L L 

) / (^ 6 0 /6 2 dx 0 ) (4.8) 

For large amplitude, the main effect of fixed wavelength is to slow the 
growth of troughs; in fact, 6 is prevented from reaching zero. In contrast, 
peak growth is only slightly accelerated. When peaks become large, then d e 
becomes negligible in comparison, as in the constant d e case, but sooner here 
since d e decreases. Thus large peaks show the same catastrophic growth to 
infinite thickness, given by (4.3a). Therefore the growth of large peaks, 
where d e is unimportant, is insensitive to (reasonable) end conditions, and 
our previous conclusion, that the growth of peaks is dramatically enhanced by 
large-amplitude effects, applies in general. 

Qualitatively this case is very similar to the previous, as shown by the 
profiles in figure 3, except that the wavelength remains constant. Figure 4 
compares the growth of the peak for fixed wavelength with that for constant 
end forces, each for an initial sinusoidal disturbance. 

4.3. LARGE- AMPLITUDE BEHAVIOR 

Where 6 >> d e , we can approximate (3.16) by 

D6/Dt = 6 2 (4.9) 

(This would apply without approximation if the layer were of finite extent, 
surrounded by fluid 2.) This integrates to: 



D6/Dt = 6 2 - 6 0 dx Q 
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<5(x 0 ,t) = 1/(1/ 6 0 (x 0 ) - t) (4.10) 

where again 6q(x 0 ) is the initial profile. Thus each cross section reaches 
infinite thickness at time t* * 3 / 6 q » this is the same growth as (4.3a). 

For this case, (3.15) simplifies to: 

6 + u x = 0 (4.11) 



so 



X 

u(x,t) = - / s(s,t) ds (4.12) 

0 



where the origin is assumed to be stationary. 

everywhere, we can use 6 dx = 6 0 dx Q , so: 

x 0 

u(x Q ) = - / 6 0 (s) ds 



As long as 6 remains finite 



(4.13) 



0 

which shows that each cross section moves at a uniform speed until the maximum 
thickness reaches the singular time. These results (4.10,4.12) apply locally 
around any large peak where d e is negligible. 



5. PLUME FORMATION 



Here we examine what happens to a thickness maximum as it approaches the 
singular time (t*). We will show how the peak becomes a plume, where the 
dense, viscous fluid drains down. We argue that the small-slope equations 
still describe the behavior of the fluid, except in an asymptotically small 
neighborhood of the plume. This is possible because that neighborhood has 
only an asymptotically small effect on the stress in the layer. Thus we can 
describe the shape of the layer around the plume, and the small-slope 
equations can be used all the way from initial conditions through the final 
draining of the layer. 

As the maximum grows, clearly the small-slope assumption must break down 
locally before the peak reaches infinite thickness; at what thickness it 
breaks down depends on the initial wavelength of the disturbance. Since the 
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initial disturbance wavelength is asymptotically long, or equivalently, the 
layer is initially asymptotically thin (e e kd 1 << 1), then when the peak 
reaches infinite thickness, that portion of the layer around the peak where 
the physical slope of the interface is 0(1) or greater will be asymptotically 
narrow compared to the whole wavelength. This follows from mass conservation, 
and from the shape near the peak approaching (as we will show) an integrable 
negative power of x (where for convenience we have chosen x = 0 at the peak). 
This point is illustrated in figure 5. So the small-slope approximation 
continues to apply almost everywhere in the layer; what is needed is a 
description of the effect of the peak or plume on the rest of the layer 

In a companion paper, we give a large-slope analysis appropriate to the 
region around a peak where the small-slope approximation no longer applies. 
There the flow is extensional, nearly vertical, driven by buoyancy and limited 
by normal viscous stresses. We find that large-slope effects do not slow down 
the growth, they only affect the detailed shape of the peak. Specifically, we 
find that the catastrophic behavior described by (4.3a) still applies (except 
for a numerical coefficient close to 1). Physically, there is nothing to 
prevent the fluid from flowing down. In this way, the peak extends to become 
a plume, where the viscous fluid continues to drain down. Of course, at some 
point the fluid will reach the lower boundary or the extending peak will 
become so long that the viscous resistance of fluid 2 becomes important. In 
the former case, only the details of the plume shape are changed, but in the 
latter, the driven flow in fluid 2 could affect the dynamics of the upper 
layer. 

To the rest of the layer the plume appears as an isolated singularity, a 
sink of fluid. The horizontal force balance must still apply even to the 
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plume, and so the reduced force in the layer F(t) is continuous across the 
plume. Recall that F(t) is an integral property of the whole layer; for a 
fixed-wavelength disturbance the dimensionless form is given by (4.7). 

Because the region where the small-slope approximation fails is integrable (as 
it must be, since no new mass is created) and asymptotically thin, its effect 
on F(t) is negligible. Equation (4.7) shows that where 6 is largest has the 
smallest effect on F. (Even after fluid starts draining out of the layer, the 
numerator in (4.7) still just gives the volume of the layer, regardless of its 
shape . ) 

As an example, consider what happens near the maximum when the initial 
disturbance is a small-amplitude sinusoid. (For simplicity, we assume the 
constant-force end conditions, but as F has little effect on a large peak the 
results apply to more general end conditions.) The previous solution (4.6g,h) 
shows that the peak becomes singular as A = be 2t /2 approaches unity. Then in 
(4.6h) the argument of the inverse tangent approaches zero (for x Q away from 
7 t ) , so near the singular time (4.6h) becomes: 

x * 2tan(x 0 /2) - x 0 (5.1a) 

for (1-A) << 1, and (1-A) 1/3 << (tt-Xq). Then near the peak (x 0 << 1): 

x * x Q 3/ 1 2 (5.1b) 

Setting A = 1 in (4.6g) gives the thickness profile at ‘the singular time: 

6 * cot 2 (x q / 2) (5.2a) 

and near the peak: 

6 * (x q /2)-2 * (3/2 x)-2/3 (5.2b) 

This shows that at the singular time the peak becomes proportional to an 
integrable negative power of x. 
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In fact, any smooth initial peak results in the same power of x. For 
x 0 << 1 , the smooth (symmetric) initial peak is asymptotically: 

6 0 « 1 + b(l - cx Q 2/ 2 ) + O(bx 0 4 ) (5.3) 

where b is the (small) initial amplitude and c is some positive constant. 

Then from (4.2a): 



6 k 1 + A ( 1 - c x n 2 / 2 ) 

1 - A(1 - c x Q 2/2) 

where again A(t) = be2t/2. Near the singular time, A + 1 and: 

6 -v 

(1-A) + c x q 2/2 



(5.4a) 



(5.4b) 



From continuity, dx/dx 0 = 6 0 / 6 <= 1/6 (for small initial amplitude), so for 
Xq << 1 and as A + 1 : 

x -V (l-A)x 0 /2 + c x 0 3/12 (5.5) 

At the singular time, A = 1, and near the singularity: 

6 - 4/(c x Q 2) » c-1/3 (3 x/2)"2/3 (5.6) 

For the sinusoidal disturbance, c = 1 and (5.6) reduces to (5.2b). 

As the singular time is approached, the peak can be described by a 
similarity solution. Defining n(x 0 ,t) by: 

(1-A(t))n 2 i c x q 2/2 (5.7a) 

reduces (5. 4b, 5. 5) to: 

f ( n) = ( 1-A) 6/ 2 = l/(l+n 2 ) (5.7b) 

£(n) = 3/(2c)(l-A)-3/2 x = n 3 + 3n (5.7c) 

This is a particular case of the general similarity solution for large 6 given 
in appendix D. The comparison is clarified by noting that: 

A = e 2 ^, T = t-t* (5.8a) 



so that as A + 1 , t + 0, then to leading order in t: 
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(1-A) «► 2t (5.8b) 

The general similarity solution shows that a peak of the more general form 
6 0 * 1 + b ( 1 - c | x 0 | n ) , where n > 0, gives a singularity of the form 
6(x,t*) a |x|~ m , where m = n/(n+l). In other words, any (reasonable) initial 
peak becomes , at t*, a singularity proportional to an integrable negative 
power of x. 

After the peak becomes singular (i.e., a plume) the shape of the layer 
around the plume can be calculated using the large-amplitude forms of the 
equations given in section 4.3. For simplicity, we assume 6 is symmetric in 
x, decreasing monotonically away from the origin, and we only consider x > 0. 
Also, we choose the time origin such that the "initial" amplitude 6 0 >> d e in 
the region of interest. Then the large-amplitude equations apply: 

6(x 0 ,t) = 1/(1/ 6 0 (x 0 ) - t) (A. 10) 

X 

u ( x , t ) = - J 6(s,t) ds (4.12) 

0 

The peak becomes infinite at t* = 1/6 q( 0). Thereafter, the layer must 
move in such a way that each fluid cross section reaches the origin at the 
same time that it reaches infinite thickness. In other words, if for t > t* 
we designate the fluid element x Q at the singularity by x Q = x*(t), so 
x*(t*) = 0, then 6o( x *(t)) = 1/t, and since <5 0 is an invertable function by 
the assumptions above, this uniquely defines x*, and thus also determines the 
strength ( 6 0 (x*)dx*/d t ) of the sink at the singularity, as a function of time. 
With these large-amplitude equations we can determine the shape around the 
singularity for all time. 

Now we consider the fate of the inverse power of x singularity that forms 
at t*r 



at t = 0: 6 0 (x 0 ) = 6(x=x 0 ,t=t*) = b x 0 "^, 0 < a < 1 



(5.9a) 



19 



where we have chosen the reference time at t*, when T = t - t* = 0, and b is 
some positive amplitude. Then the fluid cross section x c at the origin is 
specified by x Q = x*(t), where: 

x*(t) = ( bi) 1 /a (5.9b) 

Using (4.10) shows that: 

6(x 0 , t ) = b/(x Q a - x *a) (5.9c) 

so: 

x 0 

x(x 0 , T ) = / 6 0 (s)/6(s,t) ds 

x* 

= (x Q - x*) - (x*a/(l-a) ) (x 0 l"a - x*l~a) (5.9d) 

Far away from the plume (x Q » x*), the profile is still the starting profile 
from t = 0 (x « Xq, 6 ~ bx"cc) . However, very close to the plume: 

x 0 = x*(l+e), e « 1 ( 5 . 9e ) 

x * ax*e 2 /2 (5.9f) 

6/b « l/(ae:x* a ) « ( b t) 1 1 ( 2a)- 1//( 2otx ) (5.9g) 

This shows that asymptotically near the plume, the singularity goes like 
l//x, with a scale that varies in time. This asymptotic x dependence is 
actually independent of the starting conditions (as shown by Canright, 1937, 
App. C). From (5.9g) it is clear that whether 6 near the plume grows or 
shrinks is determined by whether a < 1/2 or a > 1/2, respectively. The 
special case a = 1/2 gives a steady solution. For a >1/2, the fluid drains 
away down the plume faster than it comes in from the sides, and the square- 
root singularity diminishes with time as it spreads out, to match onto the 
nearly undisturbed profile x~ a . This would be the eventual fate of an 
initially (t = 0) smooth maximum, which gives a = 2/3. Conversely, if 
a < 1/2, the square-root singularity grows as it spreads, fed from the sides 
faster than it can drain fluid away. (To get a < 1/2 would require a cusp- 



like initial maximum, which may not be physically realistic.) This solution 
(5.9a-g) is again a particular case of the general large-amplitude similarity 
solution of appendix D. To see this, define ti(xq,t): 



20 



n( x 0 , t) e x 0 /x*( t ) , n > 1 
f (n) = t 6 = 1 / ( n a - 1) 

?(n) = x/x*( T ) = n “ ( f|l” a -a)/ ( l*a) 



(5.10a) 



(5.10b) 



(5.10c) 



With the above description of how a plume first forms and how it behaves 
afterward, the small-slope equations can be used to follow the development of 
the instability from initial conditions through rapid large-amplitude growth 
all the way to the draining away of the fluid down the plumes. The results 
will be inaccurate wherever the physical slope of the interface is not small, 
but such regions comprise only a small fraction of the domain. The only 
assumption is that a plume does not exert any net horizontal force on the 
surrounding layer. This assumption may break down if a plume’s length becomes 
so much greater than the initial wavelength that the flow it drives in fluid 2 
becomes dynamically significant. 

6. CONCLUSIONS 

The central concern of this work is the nonlinear interactions between 
buoyant forces and normal viscous stresses that occur in a buoyantly unstable 
viscous layer under a shear-free horizontal boundary and over a much less 
viscous fluid, for long-wave disturbances in the range where the lower fluid 
is dynamically unimportant. After the initially uniform thickness of the 
layer has been perturbed slightly, the early growth of the (small) 
perturbation is exponential; the perturbation keeps its shape while it grows. 
But when the nonlinearities become important, the thicker parts of the layer 



thicken more rapidly while the thinner parts thin more slowly, giving in 
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general sharp peaks with broad, flat troughs in between. The accelerating 
growth of peaks leads to infinite thickness at some time t* (which depends on 
initial conditions), and the final catastrophic growth of the peak thickness 6 
is algebraic: 6 25 y/ Apg( t*-t ) . (In fact, this same sort of inverse-time 
catastrophic growth is also predicted for a power-law fluid, though the 
coefficient is different.) This shows how large-amplitude growth is 
fundamentally different from small-amplitude growth; large-amplitude effects 
dramatically enhance the growth of peaks. 

These are the results of a small-slope analysis; of course, where the 
peaks become very large the small-slope approximation breaks down and there 
the flow becomes fully two-dimensional. But as a companion work that employs 
a large-slope analysis will show, the growth of the disturbance at large 
slopes is essentially the same as that predicted here, with catastrophic 
inverse-time peak growth; although the details of the peak shape do change, 
this only changes the prediction by a numerical factor of order 1. 

Furthermore, the small-slope equations will continue to apply to the 
layer even after the formation of downwelling plumes, except in an 
asymptotically narrow neighborhood around each plume. This is possible 
because the plumes do not change the horizontal force balance (unless the flow 
they drive in the lower fluid becomes dynamically significant). Applying the 
equations up to the singular time shows that the first singularity should have 
the local shape 6 « |x|“ 2/3 , but that afterwards, as the plume drains the 
layer, the singularity changes shape to 6 « |x|~ 1/2 . This behavior is 
clarified by a family of similarity solutions, appropriate where 6 is large. 

Thus the small-slope analysis can be extended to describe the development 
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of the layer for all time, except in those narrow regions where the physical 
slope of the interface becomes large. 

APPENDICES 

A. LINEARIZED SOLUTION 

Here we calculate initial growth rates for the general problem with 
arbitrary wavelength, depths, and viscosities. Initially the perturbations 
are small, so the problem can be linearized by applying the interface 
conditions at z = dj and assuming the slope remains small, but of course 
including the buoyant pressure term due to the perturbations. (This is an 
extension of the analysis given by Whitehead and Luther, 1975, to include the 
effects of finite depth in fluid 2.) Thus the interface conditions (2.3c-f) 
become : 



at z = dj : [u] = [u(u 2 +w x )] = 0 


(A. la) 


[~P + 2yw z ] = ApgU-dj) 


(A. lb) 


St = w 


(A.lc) 



where again [] indicates the jump in value from fluid 1 to fluid 2. Then the 
equations are separable, and a simple analytical solution is possible. The 
perturbation is assumed to be sinusoidal, but as the problem is linear, any 
(small) initial condition can be constructed by superposition. 

The solution is given below, where the subscripts 2 and 1 refer to the 
two fluids, k is the wave number, a(t) is the dimensionless amplitude, 

Z e z - (d 1 +d 2 ) is the coordinate in fluid 2, a = y 2 /yi * s t * le viscosity 
ratio, and in the coefficients A, B, E, F, and their common denominator D 
these abbreviations are used : k e 2kd } , K e 2kd 2 , c = cosh(fc), C e cosh(K), 
s e sinh(k) , S e sinh(K). 
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6 = dj [1 + a ( t ) cos (kx ) ) (A.2a-j) 

"l*! = ( Apgd i / p 2 ) a(t)/2 sin(kx) [A sinh(kz) + B kz cosh(kz)] 

= (Apgdj/i^k 2 ) a(t)/2 sin(kx) [E sinh(kZ) + F kZ cosh(kZ)] 

Pi = (Apgd^ { -a(t) B cos(kx) cosh(kz) } 

P 2 = (Apgd^ { -a(t) a F cos(kx) cosh(kZ) } 

A = (-1/D) { ( 2( S-K) + ak(C-l)) sinh(kd 1 ) + 

Ik(S-K) + a (kK + 2 ( C— 1 ) ) ) cosMkdj) } 

B e (2/D) { (S-K + aK) sinh(kd!) + a (C-l) cosh(kd 1 ) } 

E = (1/D) { (2a(s-k) + K(c-l)] sinh(kd 2 ) + 

[aK(s-k) + Kk + 2 ( c— 1 ) ] cosh(kd 2 ) } 

F e (-2/D) { [a(s-k) + k] slnh(kd 2 ) + (c-1) cosh(kd 2 ) } 

D E ( S-K) (s+k) + 2 a ( Cc-l+Kk ) + a 2 (S+K)(s-k) 

Then from = w(z=d 1 ) we get the growth rate: 



a(t ) = a(0) e^t 


(A. 3a) 


o = ( Apgd 1 j ) a 


(A. 3b) 


1 ( S-K) ( c-1 ) + a( s-k) (C-1 ) 

o = — 


(A. 3c) 



k ( S-K) ( s+ic) + 2a(Cc-l+Kk) + a 2 (S+K)(s-k) 



The symmetry of the problem is apparent in the solution. Had we chosen 
fluid 2 as the reference rather than fluid 1, the solution would have the same 
form. Thus we can, without loss of generality, assume that fluid 2 is the 
deeper layer: K > k. 

This solution is governed by three independent dimensionless parameters: 
the non-dimensional wavenumbers, or depths, k and K (or equivalently k and the 
depth ratio 6 = d 2 /d 1 ), and the viscosity ratio a e y 2 /Mi. It should be noted 
that o is a monotonically decreasing function of a, i.e., if we increase p 2 
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while ui stays constant the growth rate can only decrease. Also note that if 
8 -► 00 then (A. 3c) reduces to: 

5-1 (c-1) j (A3d) 

k (s+k) + 2ac + a 2 (s-k) 

which is just the result of Whitehead and Luther (1975). Figure 6a shows the 
effects of viscosity ratio on a(k) for 8 = figure 6b shows the effects of 
finite depth for a << 1. 

When £> » 1 there are well defined regimes of growth where different 
force balances are dominant. These are shown schematically, with the 
corresponding growth rates, in figure 7. As we discuss the various growth 
mechanisms below, it should be kept in mind that the same mechanisms continue 
to apply as long as the slope of the interface remains small, which for long 
waves includes large-amplitude growth. (This point is supported by the 
scaling argument in Canright, 1987, App. B.) 

For sufficiently short waves (k >> min( 1 ,max( 8“ 1 , a* 1 / 3 ) ) ) the disturbance 
sees neither boundary and a -► [k(l+a)] -1 , so if one viscosity is much larger, 
that one limits the growth. The dominant mechanism here is that the vertical 
motion of the interface is resisted by normal viscous stresses in the more 
viscous fluid, and the growth rate diminishes with decreasing wavelength. 

At the other extreme, for sufficiently long waves 
(k << min(8“ 1 ,/( 6/a) ,/(a6) ) ) the boundaries confine the flow to be mainly 
horizontal, limited by shear stresses at the interface. Then 
a -► k 2 ( l+8/a)/12 , and the controlling viscosity depends on if a > g or not. 

Between these extremes, the waves are long compared to layer 1 so the 
motion of the interface is primarily horizontal, and there are four different 
regimes. In one (cT 1 << k << cT 1/3 , 1 « a << 6 3 ), fluid 2 is relatively 
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immobile and the growth is limited by the shear across layer 1, giving the 
same growth rate as the previous case, i.e., k 2 / 1 2 . (This possibility was 
apparently overlooked by Whitehead and Luther.) In another 
(8 _1 << k << minCajcT 1 )), the slight resistance of fluid 2 (the rate- 
controlling viscosity is y 2 ) gives a small shear gradient across layer 1, 
which over the long wavelength is sufficient to balance the buoyancy, giving 
a + k/4a. 

In the other two regimes, the less viscous fluid is effectively passive, 
and the primarily horizontal motion of the more viscous layer is limited by 
normal viscous stresses. The resulting growth rate is nearly independent of 
wavelength, and includes the maximum growth rate possible for a given 
viscosity contrast a where this behavior occurs. (For other a, i.e., 

1 << a << we expect the fastest growth at the crossover between short- and 

long-wave behavior, at k max * 2.9/a 1/3 , Om ax ~ 0.23/a 2/3 .) When fluid 2 is 
much more viscous (a >> 8 3 )> this regime (/(8/a) << k << 8~ 1 ) gives a 8/4a, 
while for fluid 1 more viscous ( a << 1) this regime (max ( a , /( a8) ) << k << 1) 
gives a + 1/4. 

In the last case, consideration of higher-order effects shows that, for 
8“ 5 << a << 1, a * ( 1 / 4 ) ( 1 - (a/k + k 4 /720)). This broad maximum peaks at 
k m ax 155 (180a) 175 = 2.8 a 175 and o max K ( 1 / 4 ) ( 1 - 0.44 d 475 ). As an indication 
of the flatness of this peak, for the range a 475 < k < 5.2 a 1720 , a > ( 1 / 4 ) ( 1 
- a 1/5 ). When a << B“ 5 , the finite depth modifies the maximum growth rate 
giving a K ( 1 / 4 ) ( 1 - (3a/8k 2 + k 4 /720)) with a broad peak at 
k max - 3.2(a/8) 176 and 5 max - ( 1/4) ( 1 - 0 . 44 ( a/ 6) 2 / 3 ) . 

The present work is only concerned with the case of a thin, viscous layer 
over a less viscous, deep layer, so a << 1 and 6 >> 1 . We further restrict 
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our consideration to the mechanism giving the fastest growth, i.e., the last 
regime considered above, where a * 1/4. The lowest-order finite-amplitude 
analysis (section 3) therefore predicts that the growth is independent of 
wavelength . 

B. POWER-LAW FLUID 

The long-wave quasi-one-dimensional analysis of section 3 is not limited 
to Newtonian rheology; any constitutive relation can be accomodated. The 
important point is that both surfaces of the layer are shear-free (in the 
wavelength range where the lower fluid is passive), so throughout the layer 
shear stresses are 0(e) smaller than normal stresses, and the latter are 
independent of z to 0(e 2 ). For a one-dimensional disturbance (3.9-10) become: 



where is the deviatoric stress tensor. 

Below, we consider the particular case of a fluid with a power-law 
rheology. In mantle convection, the oceanic lithosphere makes up most of the 
cold, stiff upper thermal boundary layer. The lithosphere behaves like a 
rigid plate rather than a viscous fluid layer, so the results for a Newtonian 
fluid may be a poor model for destabilization of the lithosphere and the 
initiation of subduction. A better model might be a fluid where the stress 
depends on some power of the strain rate. This will not accurately describe 
fracturing, but can at least incorporate weakening at high rates of strain. 

For a power-law fluid in quasi-one-dimensional flow: 



Apg6 2 /2 + 2<5t xx = F(t) 
D6/Dt = S t + u6 x = -<$u x 



(B. 1 ) 



(B. 2) 



T xx ” Pr! u xl n 



(B. 3) 



where y r has the appropriate units. Combining (B.l-3) we get a single 
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(dimensionless) Lagrangian equation describing the evolution of the thickness 
of the layer as we follow a material cross section: 

D 6 /Dt = sgn( 6 -d e (t)) 6 [ | 6 2 -d e 2 ( t ) | / 6 ] m (B.4) 

where m = 1/n, <5 = 5 /h, t e (Apgh/ 8 p r ) m t, and d e (t) is the nondimens ional ized 
current equilibrium thickness as before. For the special case of a Newtonian 
fluid (n = m = 1) this reduces to (3.10). When d e (t) is specified, (B.4) can 
be integrated numerically from the initial thickness 6 q( x q)‘ (Now drop 
tildes.) 

For an infinitely long layer with a localized disturbance (i.e., constant 
end forces: d e = 1 ) and m an odd integer: 

D 6 /Dt = ( 6 2 -l ) m / 6 m ~ 1 ( B . 5 ) 

This can be integrated by parts. For example, if m = 3: 

t( 6 ) = 1/ 4[ 6 / ( 6 2 -l ) 2 + 1 / 2 1 6 / ( 6 2 - 1 ) + 1/2 In | ( 6-1 )/ (6+1 ) | ]] (B.6) 
where t( 6 ) e t*-t is the time remaining before the thickness of the fluid 
cross section goes to °o or 0. Profiles ( 6 (x) for various t) are shown in 
figure 8 for m = 3 and m = 9, from an initial sinusoidal disturbance, using 
solutions of the above form. 

The effect of increasing m is seen to be concentration of the deformation 
into narrow regions at the centers of peaks and troughs, while elsewhere the 
profile becomes linear in x. The growth remains slow initially, then suddenly 
becomes catastrophic after a certain threshold has been reached. A trough 
reaches this threshold and necks off much sooner than a peak of equal initial 
amplitude blows up. This can be seen in figure 9, which shows the growth of 
positive and negative perturbations over time. 

When the amplitude a = 6-1 is small, the approximate solution is: 
a w a 0 / 1 1 “ (m- 1 ) ( 2 a 0 ) m ~ 1 1 ] 1 1 ) 



( B . 7 ) 
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where a Q is the initial amplitude. This slow growth lasts for a period of 
roughly At ~ 1/ f (m-1 ) ( 2a 0 ) m ~ 1 ] . For large amplitudes, the growth again 
becomes algebraic in the time remaining before the singularity is reached at 
t = t* : 

6 o > 1 : 5 « l/(m(t*-t)) (B.8a) 

6 0 < 1: 6 * m(t*-t) (B.8b) 

This catastrophic growth occurs over a time scale At ~ 1/m (dimensionally 
( 8 p r / Apgh ) m /m ) . 

As for the Newtonian fluid, a disturbed region under constant end forces 

will tend to stretch out to some extent in the x direction. To keep the 

wavelength constant, the end forces must vary in time to give: 

L 

^ [ (6 2 -d e 2(t))/6] m 6 dx = 0 (B.9) 

While we have not done a numerical calculation for this case, we speculate 
that enforcing constant wavelength does not significantly alter the 
qualitative behavior, except to inhibit the layer from necking off. (When a 
fluid cross section reaches zero thickness, locally the strain must become 
infinite; for smooth initial conditions, this is incompatible with constant 
overall strain.) Also, peak growth may be slightly enhanced. 

The main point is that when peaks get large (compared to the current 
equilibrium thickness), the large-amplitude effects still produce catastrophic 
growth of the form 6 « l/(t*-t), as for a Newtonian fluid, where t*-t is the 
time remaining before the singularity. 

C. LARGE- AMPLITUDE SIMILARITY SOLUTION 

The equations appropriate to large-amplitude disturbances admit a rich 
family of similarity solutions, which illustrate a variety of behaviors. 

While such solutions demand particular initial conditions, nonetheless these 
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solutions can be interpreted as good local approximations for situations 
arising from arbitrary initial conditions. Two cases are of particular 
interest, in light of section 5: one describes how a smooth finite peak 
evolves to an infinite singularity, the other describes how that first 
singularity changes shape as the plume evolves. 

Where 6 >> d e (in dimensionless terms), we can effectively set d e = 0 to 
get the large-amplitude equations, valid locally. In fact, the same equations 
would apply everywhere if somewhere the layer has zero thickness; for the case 
of constant end forces, minima in the layer (if initially less than d e ) reach 



zero thickness in finite time. Then: 

6 t + (u6) x =0 (C. la) 

6 + u x =0 ( C . 1 b ) 

We assume a similarity solution of the following form: 

x = a( t ) £ , 6 = A(t)f(C), u = a(t )A( t)g(C) (C.2) 

Then the system (C.l) becomes: 

(A'/A2)f - (a'/(aA))$f' + (fg)' = 0 (C.3a) 

f + g ' = 0 (C. 3b) 



Similarity then requires that A'/A 2 be a constant. If that constant is 
not zero we can normalize it by choosing the scale of A, and by a suitable 
choice of time origin: 

A ( t ) = 1/t (C.Aa) 

Then, because 6 must be non-negative, solutions where f > 0 will apply for 
t > 0, when A is decreasing, and where f < 0 will apply for t < 0, when A is 
increasing catastrophically to t = 0. 

Also, a'/(aA) must be constant, say X. (As X is insensitive to the scale 
of a(t), it cannot be normalized.) Then let: 



30 



a(t) = |t|A (C.4b) 

For t > 0, positive X gives a profile that spreads out in the x direction, 
while for negative X the profile contracts; the converse is true for t < 0. 

Combining (C.3a,b) : 

(gg')' - Kg" ~ g' = 0 (C.5) 

which integrates to: 

(g-xc)g' + (x-l)g = C (C. 6) 

where C is an arbitrary constant. By a translation, corresponding to moving 
the origin to a different part of the layer moving at a different speed, the 



constant can be removed: 

Q = £ - C/(X 2 -X), q(c) = g(5) - C/(A-1) (C.7a) 

(q-Ac)q' + (A-l)q = 0 (C.7b) 

as long as x t 0,1. Then the equation becomes separable for: 

Q(c) = cq(?) 

(Q-A)cQ' + (Q-DQ = o (C.8) 

Integrating by partial fractions and substituting back gives: 

C = q + Sgn( ?-q)D| q | k (C.9a) 

f = -1/(1 + sgn(q(^-q))kD|q| k -“l] (C.9b) 



where k = X/(X-1) and D > 0 is the second arbitrary constant of integration. 
This is the general solution; though the differential equation is nonlinear, a 
phase-plane analysis verifies that this gives all solutions (for A(t) not 
constant, X * 0,1) except the trivial solutions f = 0 and f = -1. 

The solution describes a variety of behaviors, depending primarily on X 
and on which branch of the solution is chosen. The constants C and D affect 



the x origin and scale, respectively. 
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For example, consider x > 1 ( so k > 1), C = 0, D = 1, and the particular 
branch : 

5 = g + sgn(g) | g | k (C. 10a) 

f = -1 / [ 1 + k|g| k "l] (C.lOb) 

Since f < 0, this solution only applies for t < 0; it describes the growth of 
a peak up to the singular time when 6 at the peak becomes infinite. The 
asymptotics in g reveal the behavior: 

as |g| + 0: 6 + | t | “1 [ 1 - k | t | _k | x | k ~l ] , u + -x/|t| (C.lOc) 

as |^| <$ + |x| _ l/X/k, u ->• -sgn(x) | x | 1 /k (C.lOd) 

The first (C.lOc) applies for x near the origin, as long as t is non-zero. 

The peak is of fairly general shape, but for the peak to be analytic in x, 

X must be 3/2 (k = 3). The second (C.lOd) applies far from the origin early 
on (t << 0), but applies ever nearer until at the singular time, it applies 
for all x (* 0). The asymptotic shape (C.lOd) is independent of time; as the 
peak grows it fills in the integrable negative power of x shape. For a smooth 
peak, X = 3/2 and at the singular time 6 « |x|~ 2/ 3 (for any D) . 

As another example, X > 1, C = sgn(^), D = (X-l)k/A, and the branch: 

5 = g “ sgn(g)A -1 [((A-l)|g| + l) k - 1] (C. 11a) 

f = 1/ [ ( ( X-l ) | g | + l)k"l - l] (C. lib) 

which describes a plume, symmetric in x, for t ^ 0. Asymptotically: 

as |$| -► 0: 6 + 1 / [ t 1 2 /( 2 | x | ) ] , u + -sgn(x )/ ( 2 | x | ) / 1 1 2 (C.llc) 

as |g| + 6 + ( X | x | )-!/*, u -sgn(x) (X \ x | ) 1 /^/(X-l ) (C.lld) 

At the singular time (t = 0), the profile (C.lld), which is independent of 
time, applies everywhere. In other words, this example shows what happens to 
an initial profile proportional to an integrable negative power of x; by a 
different choice of C and D, the profile (C.lld) could be made to match 
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(C.lOd) exactly. For t > 0, (C.llc) applies near the origin; regardless of 
the starting x dependence, the shape around the singularity (plume) becomes 
proportional to l//|x|. If A < 2 then the square-root singularity and the 
corresponding strength of the sink at the origin decay with time, as the layer 
drains away. (This would be true for an initially smooth profile like (C.lOc) 
with A = 3/2.) Conversely, for A > 2, the singularity grows, as fluid comes 
In from the sides faster than it can be disposed of. The special case A = 2 
gives a steady plume: 

6 = l//( 2 | x | ) , u = -sgn(x )/ ( 2 | x | ) (C.lle) 

Also, there are related solutions for A=0,0<A<1, and A = 1, which have 
the same asymptotic shapes (C.llc,d), except that, for |^| + the forms for 
u are different, and for A = 0, 6 -► e~l x l/t as \t,\ <». 

Briefly, the other behaviors governed by the similarity solution are as 
follows. For t < 0, I.e., profiles growing to the singular time, there are 
four: a minimum thickness of zero like |x|^, k > 0, locally time-independent, 
that far away levels off to approach l/|t|; a profile that approaches zero as 
x ^ -oo and 1 / 1 1 1 as x + +«>; a symmetric finite minimum flanked by plumes 
(which could be extended periodically); and a plume whose sides level off to 
approach l/|t| (rather than 0). For t > 0, where the profile starts at the 
singular time and diminishes thereafter, there is only 'one other case: a zero 
minimum flanked by plumes (which could extend periodically). The special case 
of A(t) = 1, so a(t) = e^t, gives either an isolated plume, or a zero minimum 
where |6 X | + *> surrounded by plumes (possibly periodic). All cases with 
plumes show the l//|x| shape. And solutions with the same A(t) and A but 
different C and D can be “cut and pasted" together to give other similarity 
profiles; if g and g' are continuous this will produce reasonable results. 
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Figure 1. Rayleigh-Taylor Problem. A layer of fluid 1 over a layer 
of a different fluid 2, between shear-free boundaries, where p ? < Pi . 
Interface position given by z = 6(x,t). Equilibrium (unstable; depths 
of upper and lower layers are dj and d 2 * Growth rate of instability 
depends on viscosity ratio a = U 1 /U 2 * thickness ratio 6 = d ^ / d 2 9 and 
wavenumber k (dimensionless e = kdj). 




Figure 2. Viscous Layer Profiles 6(x) for various times t, from 
initial sinusoidal perturbation of amplitude for constant end- 

force conditions (d e = 1). Each profile represents one wavelength; 
the overall strain increases with time. 




Figure 3. Viscous Layer Profiles <5(x) for various times t, from 
initial sinusoidal perturbation of amplitude 10 ”^, f or constant 
wavelength conditions (d e (t) decreases with time). Similar to Fig. 2, 
except troughs not as deep, and the peak reaches 00 slightly sooner. 
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igure 4. Peak Disturbance Growth: 6 max (t)-l (logarithmic) vs. time, 
from initial sinusoidal perturbation of amplitude 10"4. 

(solid): for constant end force 
(dashed): for constant wavelength 

(dotted): exponential (linearized) growth, for comparison 
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Figure 5. When small-slope approximation fails, it does so only in a 
region of width 0(e) . 
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Figure 6a. Linearized Growth Rate: a(k) (logarithmic scales) 

(a) Infinite depth (8+ 00 ), for various viscosity ratios a. This 
shows that when the film is much more viscous than the fluid below 
( a <<l), the maximum growth rate applies over a broad range of 
wavelengths (to lowest order). 
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Figure 6b. Linearized Growth Rate: o(ic) (logarithmic scales) 

(b) Finite depth (various 6), for a = 10 -5 . This shows how 
finite-depth effects give a large-wavelength cutoff to the maximum 
growth regime for a<<l. 




Figure 7. Rayleigh-Taylor Growth Regimes. The various (dimension- 
less) asymptotic growth rates a (from Eq. A. 3c) are shown in the 
regions of the parameter plane (wavenumber k and viscosity ratio a, 
logarithmic scales) where they apply, assuming the depth ratio B>>1 . 
(Divisions between regions, shown as lines, are actually broad areas 
across which the growth rate varies continuously.) These growth 
regimes apply while the interface slope remains small, which for long 
waves (k<<l) includes large amplitudes. 



Figure 8a. Power-Law Fluid Layer Profiles: 6(x) for various times t, 
from an initial sinusoidal perturbation of amplitude 0.02, assuming 
constant end-force conditions. 

(a) for power m = 3 

(dotted): initial, t = 0 
(dashed) : t = 1 29 
(solid): t = 153 
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Figure 8b. Power-Law Fluid Layer Profiles: 6(x) for various times t, 
from an initial sinusoidal perturbation of amplitude 0.02, assuming 
constant end-force conditions. 

(b) for power m * 9 

(dotted): initial, t = 0 
(solid): t = 8.8*109 
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Figure 9a. Power-Law Fluid Disturbance Growth: S(t) (linear scales) 
following a fluid cross section, from an initial perturbation 
amplitude of 0.02, assuming constant end forces. Growth is shown for 
both positive and negative perturbations. 

(a) for power m = 3 
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Figure 9b. Power-Law Fluid Disturbance Growth: 6(t) (linear scales) 
following a fluid cross section, from an initial perturbation 
amplitude of 0.02, assuming constant end forces. Growth is shown for 
both positive and negative perturbations. 

(b) for power m = 9 



z = 6(x,y) 




Figure 10. Force balance on a differential column of fluid (z 
upward ) . 
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